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Abstract. This paper presents an application of the TwO-Dimensional CORrelation (TODCOR) algorithm to 
multi-order spectra. The combination of many orders enables the detection and measurement of the radial veloc- 
ities of very faint companions. The technique is first applied here to the case of HD 41004, where the secondary 
is 3.68 magnitudes fainter than the primary in the V band. When applied to CORALIE spectra of this system, 
the technique measures the secondary velocities with a precision of 0.6 km s -1 and facilitates an orbital solution 
of the HD 41004 B subsystem. The orbit of HD 41004 B is nearly circular, with a companion of a 19 Mj minimum 
mass. The precision achieved for the primary is 10ms _1 , allowing the measurement of a long-term trend in the 
velocities of HD 41004 A. 
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1. Introduction 

TODCOR (TwO-Dimensional CORrelation) is a two- 
dimensional correlation technique to derive the radial ve- 
locities of both components of double- lined binary spectra 
(Zucker & Mazeh 1994). It was introduced as a gener- 
alization of the (one-dimensional) cross-correlation tech- 
nique (Simkin 1974; Tonry & Davis 1979) to deal with 
the difficulties encountered in double-lined spectra, when 
the lines of the two components cannot be easily resolved. 
Assuming the observed spectrum is a combination of two 
known spectra shifted by the radial velocities of the two 
components, TODCOR calculates the correlation of the 
observed spectrum against a combination of two templates 
with different shifts. The result is a two-dimensional corre- 
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lation function, whose peak simultaneously identifies the 
radial velocities of both the primary and the secondary. 

One of the advantages of TODCOR is its ability to 
use different templates for the primary and the secondary. 
When the primary and the secondary are of different spec- 
tral types, the simultaneous use of two different templates 
utilizes all the spectral information contained in the ob- 
served spectrum of the combined system. This feature of 
TODCOR is most important when deriving the radial 
velocities of faint secondaries (Mazeh & Zucker 1994). 
However, the ability to detect faint secondaries using 
TODCOR is still limited, depending mainly on the signal- 
to-noise ratio and the number of spectral lines. Therefore, 
Mazeh et al. (2002) and Prato et al. (2002), who searched 
for faint companions, had to apply TODCOR on IR spec- 
tra, where the flux ratio is favorable for detecting cooler 
stars. 

The modern spectrographs offer another path to en- 
hance the ability to detect faint companions. Due to the 
progress in detector technology, many of the modern spec- 
trographs produce multi-order spectra. In order to en- 
hance our ability to detect a faint companion, while max- 
imizing the precision of the measured velocities, we need 
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to combine the spectral information in all the relevant or- 
ders. Originally, TODCOR was devised to analyse only 
single-order spectra, and further generalization is there- 
fore needed in order to use the information in multiple 
orders. 

When applying TODCOR to each order separately, 
the weak signal of the secondary may produce a certain 
local peak of the correlation function at the correct sec- 
ondary velocity. However, this peak can be easily topped 
by spurious random peaks. This rules out, for example, 
calculation of the radial velocity independently for each 
order and then averaging the velocities, since many of 
the secondary velocities would have resulted from wrong 
peaks. Concatenation of the spectral orders to one single 
spectrum would require special treatment to the gaps and 
overlaps between adjacent orders. Co-adding overlapping 
regions and bridging the gaps both require interpolation, 
which would introduce artificial noise into the analysed 
spectrum. 

The approach we suggest here is to calculate the corre- 
lation function for each order separately, and then combine 
the correlation functions of all the orders. The combina- 
tion emphasizes the relevant correlation peak, and aver- 
ages out the spurious random peaks. A simple average 
of the correlation functions may not be efficient enough, 
since the combination scheme has to consider the differ- 
ent spectral information in the different orders, and weigh 
them accordingly. Zucker (2003) introduces such a scheme, 
based on a few plausible statistical assumptions. 

Using this scheme, multi-order TODCOR is applied 
in this work to the CORALIE spectra of HD 41004. As 
Santos et al. (2002) have shown, these spectra are com- 
posed of the spectra of the two visual components of the 
system, that are separated by 0.5". Santos et al. found 
that the radial velocities measured by CORALIE for these 
spectra showed a minute periodical variability. They sug- 
gested that this apparent variation was actually related 
to a much larger variation in the velocity of the faint B 
component. The presence of the spectrum of HD 41004 B 
caused a variable asymmetry in the line shape of the com- 
posite spectra, which was reflected in a minute variability 
of the measured radial velocity. 

TODCOR, when applied to the multi-order CORALIE 
spectra of HD 41004, derived the radial velocities of both 
components without any assumptions regarding the orbit. 
The derived velocities confirm the conjecture of Santos 
et al.: the B component indeed moves periodically with a 
1.3-day period, and the A component shows a long-term 
trend. We present here an orbital solution for the B com- 
ponent and derive the slope of the long-term trend of the 
A component, based on the precise velocities derived by 
TODCOR. 

The next section briefly reviews the previous results 
concerning HD 41004. The analysis and its results are pre- 
sented in Sections 3 and 4. Section 5 discusses the results 
and their implications. Section 6 concludes the paper with 
a few remarks. 



Table 1. Stellar parameters of HD 41004 A (reproduced from 
Santos et al. 2002) 



Parameter 



Value 



Spectral Type 
Parallax [mas] 
Distance [pc] 
m v 

B-V 

T off [K] 

logfl [cgs] 
M v 

Luminosity 

Mass [Mq] 

log 7?hk 

Age [Gyr] 
P rot [days] 
v sin i [km s 
[Fe/H] 



Lq] 



K1V/K2V 

23.24 

43 

8.65 

0.887 

5010 

4.42 

5.48 

0.65 

~ 0.7 

-4.66 

1.6 

~ 27 

1.22 

-0.09/+0.10 



2. Characteristics of HD 41004 

HD 41004 is a visual double system, consisting of a K1V- 
M2V pair. According to the Hipparcos catalogue, the pair 
is a common-proper-motion pair, with a V-magnitude dif- 
ference of 3.68 and a separation of 0.541" ± 0.033". The 
basic stellar parameters of HD 41004 A are summarized in 
Table 1, reproducing Table 1 of Santos et al. (2002). The 
effective temperature, surface gravity and metallicity were 
calculated by Santos et al. through Stromgrcn photome- 
try, which also showed that the star is photometrically 
stable within the instrumental precision. An independent 
estimate of the metallicity was obtained from analysis of 
the Cross-Correlation Function of CORALIE. The two in- 
dependent estimates of the metallicity are quoted in the 
table. 

The principal result of Santos et al. (2002) is the de- 
tection of a radial-velocity periodical variability in the 
CORALIE spectra of HD 41004. The variability pattern 
presented in their paper is consistent with the presence 
of a planet orbiting HD 41004 A with a 1.3-day period. 
However, Santos et al. rejected this possibility based on 
the bisector shape analysis of the CORALIE spectra. This 
analysis revealed a periodic line shape variation, having 
the same period as the radial-velocity variation. 

The interpretation Santos et al. suggested is the pres- 
ence of an object orbiting the M2V star HD 41004 Ba. In 
this model A (the K star) and B are orbiting each other 
in a wide orbit, while B is further composed of two ob- 
jects, Ba (the M star) and Bb, in a very close orbit, with 
a period of 1.3 days. According to this interpretation, the 
CORALIE spectra of HD 41004 are a combination of the 
A and B components, because their separation is much 
smaller than the diameter of the CORALIE fiber. Since 
the flux of the B component in the relevant wavelength 
range is only about 3% of the flux of the primary, large 
radial-velocity variations of B cause only small variations 
of the measured radial velocity, but their effect is mani- 
fested in the bisector shape. 
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Santos ct al. performed some simulations to test their 
interpretation and concluded that their results are con- 
sistent with the presence of a brown dwarf orbiting 
HD 41004 Ba, with a radial- velocity amplitude of the or- 
der of 5kms _1 . They also found a long-term linear trend 
in the radial velocities, which did not agree with the 
older measurements. They ascribed it to the motion of 
HD 41004 A around HD41004B, but did not rule out an 
additional component that may be involved in this varia- 
tion. 

3. Analysis 

The data we analysed comprised the 86 CORALIE spec- 
tra of HD 41004 used by Santos et al. (2002), which were 
obtained between November 2001 and February 2002, ex- 
cept for one spectrum obtained in December 2000. The 
system was further monitored and in this paper we add 
17 spectra obtained as of March 2002. 

CORALIE is a fiber-fed, cross-dispersed echelle spec- 
trograph, mounted on the 1.2-m Leonard Euler telescope 
at La Silla (Qucloz ct al. 1999). With a resolution of 
A/AA = 50 000, it covers the wavelength range 3800- 
6900 A, with 68 echelle orders. The CORALIE system uses 
a software code that produces automatically the radial ve- 
locities, but for this work we used the reduced spectral 
orders. 

The analysis used 32 CORALIE orders within the 
spectral range 4780-6820 A, after having excluded the or- 
ders that are heavily polluted by telluric lines. Bluer orders 
were excluded because the secondary signal is expected 
to be too weak in these orders. Multi-order TODCOR 
was used to combine the correlation functions from the 
32 orders. Note that the pixcl-to-radial- velocity scale is 
different for each order. Therefore, before combining the 
correlation functions we had to interpolate and re-sample 
them into a prc-dctcrmincd scale. 

TODCOR requires two spectral templates as simi- 
lar as possible to the expected primary and secondary 
spectra, and an assumed value for the flux ratio be- 
tween the two spectra. We searched for the best tem- 
plates in two datasets of spectra obtained by ELODIE 
and CORALIE. The CORALIE templates, appropriate 
for solar-type stars, were obtained as part of a program 
to derive precise abundances of planet-hosting and non- 
planet-hosting stars (Santos et al. 2000, 2001). The tem- 
plates for the later stars were obtained using ELODIE as 
part of a program studying the binarity of close M dwarfs 
(Delfossc ct al. 1998, 1999). The flux-ratio information 
was taken from Pickles (1998), using his measurements of 
typical stellar spectral energy distribution (SED), normal- 
ized according to the V-magnitude difference of HD 41004. 
By convolving the templates with a rotational broadening 
profile (e.g., Gray 1976), we created an additional degree 
of freedom which expanded our template library, allowing 
a better fit of the templates to the observed spectra. 

In all the configurations of two templates that had 
spectral types similar to K1V and M2V, the secondary ra- 



Table 2. Known stellar characteristics of the templates finally 
chosen 



Primary Secondary 
Parameter Template Template 



Name 


HD 52698 


GJ 393 


Spectral Type 


Kl 1 


M2.5 1 


B — V 


0.89 1 


1.52 1 


Tea [K] 


5235 


N/A 


log g 


4.69 


N/A 


[Fe/H] 


0.21 


N/A 
< 2.9 2 


vs'mi [kms J ] 


2.0 



1. Hiinsch et al. 1999 



2. Delfosse et al. 1998 

dial velocities exhibited a clear nearly circular orbit with 
an amplitude of around 6 km s^ 1 , and a linear trend of the 
primary. We finally chose the configuration which yielded 
the smallest residuals in the secondary velocities, relative 
to the best-fit orbital solution. Table 2 lists the known stel- 
lar characteristics of HD 52698 and GJ 393 - the two tem- 
plates finally chosen. The stellar parameters of HD 52698 
(effective temperature, surface gravity and metallicity) 
were derived as in Santos et al. (2001), whereas an es- 
timate of the v sin i was obtained using the calibration of 
the CORALIE Cross-Correlation Function width (see the 
Appendix in Santos ct al. (2002)). 

Tabic 2 quotes the original equatorial rotational ve- 
locities for both templates. Even mild broadening of the 
primary template (usinz = lkms -1 ), on top of its orig- 
inal t>sini of 2kms~ 1 , degraded the solution consider- 
ably. On the other hand, the secondary template could be 
broadened by vsini up to 5kms _1 without significantly 
affecting the solution. Finally we used a broadening of 
vsini — lkms -1 , which yielded the smallest residuals of 
the secondary velocities. 

Figure 1 shows the flux ratio as a function of wave- 
length, corresponding to the chosen spectral types of K1V 
and M2.5V, using the SED library of Pickles (1998). The 
dashed line in the figure shows the flux ratio obtained by 
assuming a black-body radiation law for the two spectra, 
with a temperature of 3500 K for the secondary. It is clear 
that while the black-body model roughly fits the detailed 
flux-ratio, it is not accurate enough and there are large dif- 
ferences, e.g., around 6500 A. In any case, using a constant 
flux ratio as in the single-order application of TODCOR 
is clearly not sufficient for such a wide spectral range. 

Figure 2 demonstrates the way multi-order TODCOR 
measures the radial velocity of HD 41004 B. The Figure 
presents one-dimensional "cuts" of the two-dimensional 
correlation function. In each "cut" the correlation is shown 
as a function of the secondary velocity, while the pri- 
mary velocity is fixed according to the location of the 
two-dimensional maximum. Because of the small flux- 
ratio between the two templates, we expect the correla- 
tion value to change very little when the secondary tem- 
plate alone is shifted. However, in each "cut" we still ex- 
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JD = 2452312.722 



JD = 2452343.590 



Wavelength [Angstroms] 



Fig. 1. The flux ratio in the relevant wavelength range, calcu- 
lated according to Pickles (1998). The dashed line represents 
the flux ratio corresponding to black-body radiation laws of 
3500 K and 5240 K. 



pect to see a relative peak around the correct secondary 
velocity. The Figure is divided into two columns, corre- 
sponding to two different exposures (JD = 2 452 312.722 
and JD = 2 452 343.590). Using the best-fit orbit, the ex- 
pected secondary velocity for the spectrum used in the left 
column of the Figure is 35. 2 km s -1 , while for the right 
column it is 42.2 km s _1 . These velocities are marked on 
the figure by dashed lines. The primary velocity is repre- 
sented by a dotted line in the left column, while in the 
right column it almost coincides with the secondary ve- 
locity. For this graphic demonstration only, a third-order 
best-fit polynomial was subtracted from all the plotted 
functions in order to accentuate the local peak. 

The three upper panels in each column demonstrate 
the problems in the single-order TODCOR. On the left, 
the first upper panel exhibits a very prominent peak in a 
wrong velocity. The second and third panels show a mod- 
erate peak at about the expected velocity In almost all 
32 orders, some local peaks appeared around the correct 
velocity. The lower panel shows the result of combining 
the correlation functions of all the analysed orders. The 
correct peak is clearly emphasized relative to the spurious 
peaks. 

Comparing the two columns of Figure 2, we see that we 
cannot know in advance which orders present the correct 
peak. Thus, we have to combine all 32 orders in order to 
have the correct peak emphasized. The right column of the 
Figure demonstrates another advantage of TODCOR: in 
this spectrum the secondary velocity is almost identical to 
the primary velocity. In the conventional one-dimensional 
cross-correlation, there is no way to measure the two ve- 
locities, due to the blending of the two correlation peaks. 
TODCOR, by using two different templates, allows the 
measurement of both velocities. 







Fig. 2. The upper three panels in each column show "cuts" 
(see text) of the two-dimensional correlation function for three 
selected orders. The dashed lines represent the best-fit sec- 
ondary velocities while the dotted lines represent the primary 
velocity. It is not shown on the right column where it almost 
coincides with the secondary velocity. The lower panel in each 
column shows the "cut" of the function obtained after combin- 
ing the corresponding correlation functions of all the 32 orders. 
A third-order best-fit polynomial was subtracted from all the 
shown functions to accentuate the local peak. 



The listing of the radial velocities of A and B and the 
corresponding times, can be obtained at the CDS. 



4. Results 

Figure 3 shows the resulting orbit of HD41004B, while 
the orbital elements are summarized in Table 3. Applying 
multi-order TODCOR to many orders (32), combined 
with the large number of measurements (103) yielded a 
very precise orbital solution. Thus, the radial-velocity am- 
plitude (K) was found with a precision of about 1%. This 
fine precision allowed also a very accurate estimate of the 
orbital eccentricity (e). Although very small, the eccentric- 
ity is still non- vanishing, with a significance level of 2-10 
according to the Lucy & Sweeney (1971) test. Assuming a 
mass of ~ 0.4 Mq for HD41004Ba, the companion min- 
imum mass is 19 Mj. The uncertainty of 0.25 Mj does 
not take into account the uncertainty in the mass of 
HD 41004 Ba - M Ba - A 20% uncertainty in M Ba would 
result in a 2.3 Mj uncertainty in M B b,min- 

Figure 4 shows the velocities of HD 41004 A. The ve- 
locities measured after JD = 2,452,200 show a clear 
linear trend. The best-fit line has a slope of +105 ± 
10ms" 1 year -1 , and is also shown in the Figure. As Santos 
et al. (2002) have already noticed, the first isolated mea- 
surement does not agree with the linear trend implied by 
the other measurements. A Lomb-Scargle pcriodogram of 
the de-trended velocities (Figure 5) shows no hint of the 
1.3-day periodicity, again proving that only the B compo- 
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Fig. 3. Phased radial velocities of HD 41004 B 



Table 3. Best-fit orbital solution of HD 41004 B 



p 


[days] 


1.328199 


± 


0.000081 


T 


[JD] 


2 452 339.212 ± 


0.040 


e 




0.065 


± 


0.014 


7 


[kms^ 1 ] 


41.631 


± 


0.057 




[°] 


171 


± 


11 


K 


[kms^ 1 ] 


6.192 


± 


0.081 


asini 


[KT 3 AU] 
[1O- 6 M ] 


0.7544 


± 


0.0098 


f(m) 


32.5 


± 


1.3 




[Mj] 


18.64 


± 


0.26 


N 






103 




CO-C 


[kms" 1 ] 




0.56 





42.62 - 
42.6 - 
^42.58 - 
'to 42.56 - 
_| 42.54 - 
^42.52 - 
42.5- 
42.48 - 
42.46 - „ - ' 



1850 1900 1950 2000 2050 2100 2150 2200 2250 2300 2350 2400 

JD-2450000 

Fig. 4. Radial velocities of HD 41004 A. The dashed line is the 
best fit to the velocities, ignoring the first isolated velocity. 





0.1 0.2 0.3 0.4 



0.6 0.7 0.8 0.9 



t Assuming M Ba = 0.4 Mq. 



frequency [day 1 ] 

Fig. 5. Lomb-Scargle periodogram of the velocities of 
HD 41004 A after subtracting the linear trend. The dashed line 
marks the frequency corresponding to the 1.3-day period of 
HD 41004 B. 



nent participates in the periodic motion, as suggested by 
Santos et al. (2002). In light of the trend in the veloci- 
ties of component A, we tried also to fit the B velocities 
with an additional trend but the resulting trend was not 
statistically significant. 

The mean radial velocity of A during the linear part is 
42.5768 ± 0.0009 km s _1 , which is very close to the centcr- 
of-mass velocity of B - the difference is only ARV = 
0.95 ± 0.06 km s -1 . However, this difference depends also 
on the estimated velocities we used for the templates, 
and therefore we adopt a conservative error estimate of 
0.1 km s" 1 for ARV. 



5. Discussion 

The radial velocities derived in the previous section sup- 
port the triple-system model suggested by Santos et al. 
for HD 41004. Our results allow, nevertheless, a somewhat 
more detailed study of the system. 



5.1. The long-term trend of the radial velocity of 
HD 41004 A 

As we have demonstrated in the previous section (see 
Figure 4), the radial velocities of HD 41004 A show a clear 
long-term variation over time. For the period November 
2001 to February 2002 this variation can be nicely fitted 
by a linear increase. Note, however, that one single point, 
at JD = 2 451 902.774, which was observed one year before 
all the others - on December 2000, deviates substantially 
from this linear fit. This point may indicate that the ex- 
trapolation of the linear approximation is no longer valid 
for the time of that point. Thus, we may actually be seeing 
the curved part of a long-term orbit, caused by a fourth 
component in the system. Such an object would have an 
orbital motion with a much shorter period than that of the 
AB system, which is at least 100 years, as inferred from 
the observed separation. In fact, follow-up measurements 
we are currently performing hint that this may indeed be 
the case, but a detailed solution of the orbit would still be 
premature. 
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5.2. The close orbit and the rotation of HD 41004 Ba 

As mentioned in Section 3, the process by which we chose 
the templates allows a very crude estimate of the equa- 
torial rotational velocity. The best-fit orbit was attained 
when the M2 template was convolved with a rotation pro- 
file of vsini — lkms -1 . Our procedure is not accurate 
enough to use this value as a true measurement of vsini. 
However, trying large values for vsini degraded the qual- 
ity of the orbital solution. For example, while the qual- 
ity of the solution remained almost unchanged for v sin i 
as large as 5kms _1 , the O — C RMS was doubled at 
vsini = lOkms -1 . Thus, we believe that the equatorial 
velocity of Ba is of the order of a few kms -1 , and proba- 
bly less than 10 kms -1 . Even this crude estimate is suf- 
ficient to conclude that v sin i is substantially lower than 
the expected equatorial velocity of ~ 20kms~ 1 , assuming 
synchronization (Santos et al. 2002). 

At such a close orbit, with a period of 1.33 days, it 
is usually assumed that synchronization and alignment 
of the spin with the orbital motion had been established 
(Santos et al. 2002). An obvious explanation of the ap- 
parently small v sin i would be small orbital inclination. 
Assuming alignment, small inclination implies that the 
mass of the unseen companion is considerably larger than 
its derived minimum value. 

5.3. The small eccentricity of the close orbit 

A binary or a star-planet system with a period as short 
as 1.33 days is naively believed to have been circularized. 
However, our radial-velocity solution shows that the close 
orbit has a finite non-zero eccentricity - e = 0.065 ±0.014, 
which has to be explained. 

Let us first estimate the timescale of circularization 
due to processes occurring in Ba. The relevant process 
is dissipation of the equilibrium tide through interaction 
with the convective envelope (Zahn 1989). We follow Rasio 
et al. (1996) and write: 

- Tc M 1 , a 

Tdrc ~ / M cnvq (l + q ) ( R> ' 

where the parameter t c « (MR 2 / L)i is the eddy turnover 
timescale, and its numerical value is about 0.5 yr. The nu- 
merical value of / is obtained by integrating the viscous 
dissipation of the tidal energy through the convective zone 
and is of order unity (Zahn 1977). In our case the tidal 
pumping period - P/2 - is much smaller than r c , and 
therefore only convective eddies with turn-over time less 
than P/2 contribute to the dissipation. The value of / 
is then reduced to (P/(2r c )) Q (Zahn 1989). The correct 
value of a is debated but generally assumed to be either 1 
(Zahn 1992) or 2 (Goldreich & Keeley 1977; Goldman & 
Mazeh 1991). We assume that the mass of the convective 
envelope - M env - is very close to the stellar mass - M 
(e.g., Chabrier & Baraffc 2000). The resulting timescale is 
about 10 10 yr, assuming a = 1, and certainly much larger 
if a = 2. At the age of 1.5 Gyr, it therefore seems that 



circularization through dissipation in Ba has not taken 
place. 

Now let us examine the possible involvement of Bb 
in circularization processes. If the object is substellar or 
maybe even a very late M-dwarf, radiative zones in the 
atmosphere may form either due to the internal physics 
(Burrows et al. 1997) or due to the external heating by the 
primary (Guillot et al. 1996). These would necessitate an- 
other approach to calculate the circularization timescale, 
probably based on dissipation of the dynamical tide (Zahn 
1977). The theory of orbital evolution through tidal dissi- 
pation is still debated, specially the mechanisms for dissi- 
pation of the dynamical tide (e.g., Claret et al. 1995), and 
therefore it is not clear whether it contributes significantly 
in the case of HD 41004 B. 

6. Conclusion 

We have presented in this work the first application of 
multi-order TODCOR to echelle spectra. The case of 
HD 41004 demonstrates the unique capabilities of this 
technique. It utilizes all the available prior knowledge 
regarding both spectral components (the different tem- 
plates), it is not limited to a fixed flux ratio (we used the 
SED to calculate it) and it incorporates optimally the data 
from all the relevant spectral orders. Based on the best-fit 
solutions of the present case, we estimate the precision of 
the secondary velocity of HD 41004 to be 0.56 km s -1 and 
that of the primary velocities to be 10ms -1 . The radial 
velocities yielded accurate orbital elements of the unseen 
companion of HD 41004 B, and an accurate measurement 
of the radial acceleration of HD 41004 A. The latter sug- 
gested an additional component may be present in the 
system. Follow-up observations, currently underway, tend 
to confirm this hypothesis. 

The combination of multi-order TODCOR together 
with the high signal-to-noise and high resolution of the 
CORALIE spectra render this analysis a very promising 
path toward expanding the database of spectroscopic bi- 
naries and multiple systems (like HD 41004). It may also 
facilitate the detection of planets in binary stellar systems, 
which are lately the focus of an increasing interest. 
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